      MODULE GWFSUBMODULE
        INTEGER,SAVE,POINTER ::IIBSCB,ITMIN,NNDB,NDB,NMZ,NN,ND2,IDSAVE
        REAL,   SAVE,POINTER ::AC1,AC2
        LOGICAL,SAVE,POINTER ::NDF,NNDF
        INTEGER,SAVE, DIMENSION(:),   POINTER ::ISBOCF,ISBOCU
        LOGICAL,SAVE, DIMENSION(:,:), POINTER ::OCFLGS
        LOGICAL,SAVE, DIMENSION(:),   POINTER ::OCLAY
        INTEGER,SAVE, DIMENSION(:),   POINTER ::ILSYS
        INTEGER,SAVE, DIMENSION(:),   POINTER ::NTSSUM
        INTEGER,SAVE, DIMENSION(:),   POINTER ::LN
        INTEGER,SAVE, DIMENSION(:),   POINTER ::LDN
        INTEGER,SAVE, DIMENSION(:),   POINTER ::NZ
        REAL,   SAVE, DIMENSION(:),   POINTER ::RNB
        REAL,   SAVE, DIMENSION(:),   POINTER ::DH
        REAL,   SAVE, DIMENSION(:),   POINTER ::DHP
        REAL,   SAVE, DIMENSION(:),   POINTER ::DHC
        REAL,   SAVE, DIMENSION(:),   POINTER ::DZ
        REAL,   SAVE, DIMENSION(:),   POINTER ::HC
        REAL,   SAVE, DIMENSION(:),   POINTER ::SCE
        REAL,   SAVE, DIMENSION(:),   POINTER ::SCV
        REAL,   SAVE, DIMENSION(:),   POINTER ::DCOM
        REAL,   SAVE, DIMENSION(:),   POINTER ::A1
        REAL,   SAVE, DIMENSION(:),   POINTER ::A2
        REAL,   SAVE, DIMENSION(:),   POINTER ::BB
        REAL,   SAVE, DIMENSION(:),   POINTER ::SUB
        REAL,   SAVE, DIMENSION(:,:), POINTER ::DP
        REAL,   SAVE, DIMENSION(:,:), POINTER ::DVB
      TYPE GWFSUBTYPE
        INTEGER, POINTER  ::IIBSCB,ITMIN,NNDB,NDB,NMZ,NN,ND2,IDSAVE
        REAL,    POINTER  ::AC1,AC2
        LOGICAL, POINTER  ::NDF,NNDF
        INTEGER, DIMENSION(:),   POINTER ::ISBOCF,ISBOCU
        LOGICAL, DIMENSION(:,:), POINTER ::OCFLGS
        LOGICAL, DIMENSION(:),   POINTER ::OCLAY
        INTEGER, DIMENSION(:),   POINTER ::ILSYS
        INTEGER, DIMENSION(:),   POINTER ::NTSSUM
        INTEGER, DIMENSION(:),   POINTER ::LN
        INTEGER, DIMENSION(:),   POINTER ::LDN
        INTEGER, DIMENSION(:),   POINTER ::NZ
        REAL,    DIMENSION(:),   POINTER ::RNB
        REAL,    DIMENSION(:),   POINTER ::DH
        REAL,    DIMENSION(:),   POINTER ::DHP
        REAL,    DIMENSION(:),   POINTER ::DHC
        REAL,    DIMENSION(:),   POINTER ::DZ
        REAL,    DIMENSION(:),   POINTER ::HC
        REAL,    DIMENSION(:),   POINTER ::SCE
        REAL,    DIMENSION(:),   POINTER ::SCV
        REAL,    DIMENSION(:),   POINTER ::DCOM
        REAL,    DIMENSION(:),   POINTER ::A1
        REAL,    DIMENSION(:),   POINTER ::A2
        REAL,    DIMENSION(:),   POINTER ::BB
        REAL,    DIMENSION(:),   POINTER ::SUB
        REAL,    DIMENSION(:,:), POINTER ::DP
        REAL,    DIMENSION(:,:), POINTER ::DVB
      END TYPE
      TYPE(GWFSUBTYPE), SAVE  ::GWFSUBDAT(10)

      END MODULE GWFSUBMODULE



      SUBROUTINE GWF2SUB7AR(IN,IGRID)
C     ******************************************************************
C     ALLOCATE ARRAY STORAGE FOR SUBSIDENCE PACKAGE.
C     READ SUBSIDENCE PACKAGE DATA.
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,      ONLY:IOUT,NCOL,NROW,NLAY,ISSFLG,NPER,NSTP,HNEW,
     1                      DELR,DELC,BUFF
      USE GWFSUBMODULE,ONLY:IIBSCB,ITMIN,NNDB,NDB,NMZ,NN,ND2,IDSAVE,
     1                      AC1,AC2,NDF,NNDF,ISBOCF,ISBOCU,
     2                      OCFLGS,OCLAY,ILSYS,NTSSUM,LN,LDN,NZ,RNB,
     3                      DH,DHP,DHC,DZ,HC,SCE,SCV,DCOM,A1,A2,BB,
     4                      SUB,DP,DVB
C
      DIMENSION IFL(13)
      CHARACTER*24 ANAME(10)
      CHARACTER*200 LINE
      DATA ANAME(1) /'   PRECONSOLIDATION HEAD'/
      DATA ANAME(2) /'ELASTIC INTERBED STORAGE'/
      DATA ANAME(3) /' VIRGIN INTERBED STORAGE'/
      DATA ANAME(4) /'     STARTING COMPACTION'/
      DATA ANAME(5) /'     DELAY STARTING HEAD'/
      DATA ANAME(6) /'   DELAY PRECOLSOL. HEAD'/
      DATA ANAME(7) /'DELAY INITIAL COMPACTION'/
      DATA ANAME(8) /'DELAY INTERBED THICKNESS'/
      DATA ANAME(9) /'   MATERIAL ZONE INDICES'/
      DATA ANAME(10)/'NUMBER OF BEDS IN SYSTEM'/
      DIMENSION IBUFF(NCOL,NROW)
C     ------------------------------------------------------------------
      ALLOCATE (IIBSCB,ITMIN,NNDB,NDB,NMZ,NN,ND2,IDSAVE)
      ALLOCATE (AC1,AC2)
      ALLOCATE (NDF,NNDF)
      ALLOCATE (ISBOCF(6),ISBOCU(6))
      ZERO=0.0
C
C1------IDENTIFY PACKAGE.
      WRITE(IOUT,1)IN
    1 FORMAT(/,'SUB7 -- SUBSIDENCE PACKAGE, VERSION 7,',
     1     ' 03/31/2006',' INPUT READ FROM UNIT',I3)
C
C2------CHECK TO SEE THAT SUBSIDENCE OPTION IS APPROPRIATE
C2------IF INAPPROPRIATE PRINT A MESSAGE & STOP THE SIMULATION.
C2------ALSO, SUM TO GET THE TOTAL NUMBER OF TIME STEPS IN THE
C2------SIMULATION.
C
      NSTPT=0
      DO 12 NS=1,NPER
      NSTPT=NSTPT+NSTP(NS)
      IF(ISSFLG(NS).NE.0.AND.NS.GT.1) THEN
       WRITE(IOUT,10)
   10  FORMAT(1X,'SUBSIDENCE CANNOT BE USED IN SIMULATIONS',
     1  ' IN WHICH STRESS PERIODS OTHER THAN THE ',/,1X,
     2  ' FIRST ARE STEADY-STATE. SIMULATION ABORTED.')
       CALL USTOP(' ')
      ENDIF
 12   CONTINUE
C
C3------ALLOCATE SPACE FOR ARRAY NTSSUM, WHICH WILL CONTAIN THE TOTAL
C3------NUMBER OF TIME STEPS PRIOR TO THE CURRENT TIME STEP.
      ALLOCATE(NTSSUM(NPER))
C
C4------READ FLAG FOR STORING CELL-BY-CELL STORAGE CHANGES AND
C4------FLAG FOR PRINTING AND STORING COMPACTION, SUBSIDENCE, AND
C4------CRITICAL HEAD ARRAYS.
      CALL URDCOM(IN,IOUT,LINE)
      LLOC=1
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IIBSCB,R,IOUT,IN)
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISUBOC,R,IOUT,IN)
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,NNDB,R,IOUT,IN)
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,NDB,R,IOUT,IN)
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,NMZ,R,IOUT,IN)
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,NN,R,IOUT,IN)
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,3,I,AC1,IOUT,IN)
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,3,I,AC2,IOUT,IN)
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ITMIN,R,IOUT,IN)
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IDSAVE,R,IOUT,IN)
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IDREST,R,IOUT,IN)
      IF(AC2.EQ.ZERO) AC2=1.0
      NDF=.TRUE.
      NNDF=.TRUE.
      IF(NNDB.LT.1) THEN
       NNDF=.FALSE.
       NNDB=0
      ENDIF
      IF(NDB.LT.1) THEN
       NDF=.FALSE.
       NDB=0
       NMZ=0
       NN=0
      ENDIF
      WRITE(IOUT,50) NNDB,NDB,NMZ,NN
   50 FORMAT(/,'         NUMBER OF SYSTEMS OF NO-DELAYED INTERBEDS:',
     1 I3,/,'              NUMBER OF SYSTEMS OF DELAY INTERBEDS:',
     2 I3,/,'                          NUMBER OF MATERIAL ZONES:',
     3 I3,/,'                    NUMBER OF NODES IN EACH STRING:',I3)
      IF(IDSAVE.GT.0) THEN
       WRITE(IOUT,52) IDSAVE
   52 FORMAT(' RESTART INFORMATION WILL BE SAVED ON UNIT ', I5,
     1 ' FOR DELAY INTERBEDS')
      ELSE
       WRITE(IOUT,53)
   53 FORMAT(' RESTART INFORMATION WILL NOT BE SAVED FOR DELAY',
     1 ' INTERBEDS')
      ENDIF
      IF(IDREST.GT.0) THEN
       WRITE(IOUT,54) IDREST
   54 FORMAT(' RESTART INFORMATION WILL BE READ FROM UNIT ', I5,
     1 ' FOR DELAY INTERBEDS')
      ELSE
       WRITE(IOUT,55)
   55 FORMAT(' RESTART INFORMATION WILL NOT BE READ FOR DELAY',
     1 ' INTERBEDS')
      ENDIF
C
C4A-----ABORT IF NO LAYERS ARE SPECIFIED FOR INTERBED STORAGE
      IF(.NOT.NNDF.AND..NOT.NDF) THEN
       WRITE(IOUT,60)
   60  FORMAT(1X,'NO LAYERS WITH INTERBED STORAGE OF EITHER TYPE ',
     1  'WERE SPECIFIED IN INPUT.',/,1X,'SIMULATION ABORTED.')
       CALL USTOP(' ')
      ENDIF
C4B-----ABORT IF NO PROPERTY ZONES ARE SPECIFIED
      IF(NDF.AND.NMZ.LT.1) THEN
         WRITE(IOUT,*) ' STOPPING-- At least one property zone must ',
     &                 'be specified for delay beds.'
         CALL USTOP(' ')
      ENDIF
C4C-----ABORT IF NOT ENOUGH NODES ARE SPECIFIED
      IF(NDF.AND.NN.LT.2) THEN
         WRITE(IOUT,*) ' STOPPING-- Number of nodes in strings for ',
     &                 'delay beds (NN) should be at least 2.'
         CALL USTOP(' ')
      ENDIF
C
C5------IF CELL-BY-CELL TERMS TO BE SAVED THEN PRINT UNIT NUMBER.
   70 IF(IIBSCB.GT.0) WRITE(IOUT,80) IIBSCB
   80 FORMAT(1X,'CELL-BY-CELL FLOW TERMS WILL BE SAVED ON UNIT',I3)
C
C5A-----IF OUTPUT CONTROL FOR PRINTING ARRAYS IS SELECTED PRINT MESSAGE.
      IF(ISUBOC.GT.0) WRITE(IOUT,90)
   90 FORMAT(1X,'OUTPUT CONTROL RECORDS FOR SUB PACKAGE WILL BE ',
     1 'READ EACH TIME STEP.')
C
C6------READ IN MODEL LAYER NUMBERS FOR EACH SYSTEM OF INTERBEDS,
C6------FOR LAYERS WITHOUT DELAY.
      IF(NNDF) THEN
       ALLOCATE(LN(NNDB))
       WRITE(IOUT,100) NNDB
  100  FORMAT(/,' MODEL LAYER ASSIGNMENTS FOR EACH OF',I3,' NO-DELAY',
     1  ' SYSTEMS OF INTERBEDS:')
       CALL URDCOM(IN,IOUT,LINE)
       READ(LINE,*) (LN(N),N=1,NNDB)
       WRITE(IOUT,115) (LN(N),N=1,NNDB)
  115  FORMAT(1X,25I4)
       DO 120 N=1,NNDB
       IF(LN(N).GE.1.AND.LN(N).LE.NLAY) GO TO 120
       WRITE(IOUT,118)
  118  FORMAT(/,' IMPROPER LAYER ASSIGNMENT FOR NO-DELAY SYSTEM OF ',
     1  'INTERBEDS.',/,' ABORTING...')
       CALL USTOP(' ')
  120  CONTINUE
      ELSE
       ALLOCATE(LN(1))
      ENDIF
C
C7------READ IN MODEL LAYER NUMBERS FOR EACH SYSTEM OF INTERBEDS,
C7------FOR LAYERS WITH DELAY.
      IF(NDF) THEN
       ALLOCATE(LDN(NDB))
       WRITE(IOUT,135) NDB
  135  FORMAT(/,' MODEL LAYER ASSIGNMENTS FOR EACH OF',I3,' DELAY',
     1  ' SYSTEMS OF INTERBEDS:')
       CALL URDCOM(IN,IOUT,LINE)
       READ(LINE,*) (LDN(N),N=1,NDB)
       WRITE(IOUT,115) (LDN(N),N=1,NDB)
       DO 140 N=1,NDB
       IF(LDN(N).GE.1.AND.LDN(N).LE.NLAY) GO TO 140
       WRITE(IOUT,138)
  138  FORMAT(/,' IMPROPER LAYER ASSIGNMENT FOR DELAY SYSTEM OF ',
     1  'INTERBEDS.',/,' ABORTING...')
       CALL USTOP(' ')
  140  CONTINUE
      ELSE
       ALLOCATE(LDN(1))
      ENDIF
C
C8------ALLOCATE SPACE FOR THE ARRAYS HC, SCE, SCV, AND SUB.
      NCR=NROW*NCOL
      NND1=NCR*NNDB
      ND1=NCR*NDB
      ND2=0
C
C9-----READ IN ARRAY RNB TO SEE HOW MANY STRINGS OF NN CELLS ARE NEEDED.
      IF(NDF) THEN
       ALLOCATE(RNB(ND1))
       NNSUM=0
       DO 190 KQ=1,NDB
       LOC1 = 1+(KQ-1)*NCR
       LAYNUM=LDN(KQ)
       WRITE(IOUT,144) KQ
 144   FORMAT(/,1X,' SYSTEM',I4,' OF DELAY BEDS:')
C       CALL U2DREL(RNB(LOC1),ANAME(10),NROW,NCOL,LAYNUM,IN,IOUT)
       CALL U2DREL(BUFF(:,:,1),ANAME(10),NROW,NCOL,LAYNUM,IN,IOUT)
       CALL GWF2SUB72D1D(BUFF(:,:,1),NCOL,NROW,RNB,ND1,LOC1)
       DO 180 N=1,NCR
       IF(RNB(LOC1+N-1).GE.1.0) NNSUM=NNSUM+1
  180  CONTINUE
  190  CONTINUE
       ND2=NN*NNSUM
      ELSE
       ALLOCATE(RNB(1))
      ENDIF
      IF(ND2.LT.1.AND.NDF) THEN
         WRITE(IOUT,*) ' STOPPING-- Delay beds were not found in ',
     &                 ' array specifying numbers of delay beds (RNB).'
         CALL USTOP(' ')
      ENDIF
C
C10-----ALLOCATE MEMORY.
      ALLOCATE(OCFLGS(13,NSTPT))
      ALLOCATE(OCLAY(NLAY))
      IF(NNDF) THEN
         ALLOCATE(HC(NND1))
         ALLOCATE(SCE(NND1))
         ALLOCATE(SCV(NND1))
         ALLOCATE(SUB(NND1))
         ALLOCATE(ILSYS(NNDB))
      ELSE
         ALLOCATE(HC(1))
         ALLOCATE(SCE(1))
         ALLOCATE(SCV(1))
         ALLOCATE(SUB(1))
         ALLOCATE(ILSYS(1))
      ENDIF
      IF(NDF) THEN
         ALLOCATE(NZ(ND1))
         ALLOCATE(DZ(ND1))
         ALLOCATE(DCOM(ND1))
         ALLOCATE(DHP(ND2))
         ALLOCATE(DH(ND2))
         ALLOCATE(DHC(ND2))
         ALLOCATE(DP(NMZ,3))
         ALLOCATE(DVB(NDB,4))
         ALLOCATE(A1(NN))
         ALLOCATE(A2(NN))
         ALLOCATE(BB(NN))
      ELSE
         ALLOCATE(NZ(1))
         ALLOCATE(DZ(1))
         ALLOCATE(DCOM(1))
         ALLOCATE(DHP(1))
         ALLOCATE(DH(1))
         ALLOCATE(DHC(1))
         ALLOCATE(DP(1,1))
         ALLOCATE(DVB(1,1))
         ALLOCATE(A1(1))
         ALLOCATE(A2(1))
         ALLOCATE(BB(1))
      ENDIF
C
C11-----READ ARRAYS.
      NCR=NROW*NCOL
      ANNI=0.5/(FLOAT(NN)-.5)
C
C12-----READ RESTART RECORDS IF THIS SIMULATION CONTINUES FROM A
C12-----PREVIOUS SIMULATION
      IF(NDF) THEN
       IF(IDREST.GT.0) THEN
        READ(IDREST) NND2
        IF(NND2.EQ.ND2) THEN
         WRITE(IOUT,242)
  242    FORMAT(' HEAD AND PRECONSOLIDATION HEAD FOR DELAY BEDS ARE',
     1   ' BEING READ FROM RESTART RECORDS')
         READ(IDREST) (DH(N),N=1,ND2)
         READ(IDREST) (DHC(N),N=1,ND2)
         DO 250 N2=1,ND2
         DHP(N2)=DH(N2)
  250    CONTINUE
        ELSE
         WRITE(IOUT,252)
  252    FORMAT(' HEAD AND PRECONSOLIDATION HEAD FOR DELAY BEDS ',
     1   'CANNOT BE READ FROM RESTART RECORDS',/,
     2   ' SIMULATION ABORTING')
         CALL USTOP(' ')
        ENDIF
       ENDIF
      ENDIF
C
C13-----READ IN ARRAYS FOR SYSTEMS OF NO-DELAY INTERBEDS.
      IF(NNDF) THEN
       DO 260 KQ=1,NNDB
       K=LN(KQ)
       LOC1=1+(KQ-1)*NCR
       WRITE(IOUT,256) KQ
  256  FORMAT(/,1X,' SYSTEM',I4,' OF NO-DELAY BEDS:')
C       CALL U2DREL(HC(LOC1),ANAME(1),NROW,NCOL,K,IN,IOUT)
       CALL U2DREL(BUFF(:,:,1),ANAME(1),NROW,NCOL,K,IN,IOUT)
       CALL GWF2SUB72D1D(BUFF(:,:,1),NCOL,NROW,HC,NND1,LOC1)
       WRITE(IOUT,256) KQ
C       CALL U2DREL(SCE(LOC1),ANAME(2),NROW,NCOL,K,IN,IOUT)
       CALL U2DREL(BUFF(:,:,1),ANAME(2),NROW,NCOL,K,IN,IOUT)
       CALL GWF2SUB72D1D(BUFF(:,:,1),NCOL,NROW,SCE,NND1,LOC1)
       WRITE(IOUT,256) KQ
C       CALL U2DREL(SCV(LOC1),ANAME(3),NROW,NCOL,K,IN,IOUT)
       CALL U2DREL(BUFF(:,:,1),ANAME(3),NROW,NCOL,K,IN,IOUT)
       CALL GWF2SUB72D1D(BUFF(:,:,1),NCOL,NROW,SCV,NND1,LOC1)
       WRITE(IOUT,256) KQ
C       CALL U2DREL(SUB(LOC1),ANAME(4),NROW,NCOL,K,IN,IOUT)
       CALL U2DREL(BUFF(:,:,1),ANAME(4),NROW,NCOL,K,IN,IOUT)
       CALL GWF2SUB72D1D(BUFF(:,:,1),NCOL,NROW,SUB,NND1,LOC1)
  260  CONTINUE
C
C14-----INITIALIZE ARRAYS FOR SYSTEMS OF NO-DELAY INTERBEDS.
       DO 280 KQ=1,NNDB
       K=LN(KQ)
       NQ=(KQ-1)*NCR
C       NK=(K-1)*NCR
       DO 270 IR=1,NROW
       NQR=NQ+(IR-1)*NCOL
C       NKR=NK+(IR-1)*NCOL
       DO 270 IC=1,NCOL
       LOC2=NQR+IC
C       LOC2H=NKR+IC
C
C15------MULTIPLY STORAGE BY AREA TO GET STORAGE CAPACITY.
       AREA=DELR(IC)*DELC(IR)
       SCE(LOC2)=SCE(LOC2)*AREA
       SCV(LOC2)=SCV(LOC2)*AREA
C
C16-----MAKE SURE THAT PRECONSOLIDATION HEAD VALUES
C16-----ARE CONSISTANT WITH STARTING HEADS.
       IF(HC(LOC2).GT.HNEW(IC,IR,K)) HC(LOC2)=HNEW(IC,IR,K)
  270  CONTINUE
  280  CONTINUE
      ENDIF
      IF(NDF) THEN
C
C17-----READ IN TABLE OF MATERIAL PROPERTIES: K, Sse, Ssv FOR EACH
C17-----OF NMZ ZONES.
       WRITE(IOUT,295)
  295 FORMAT(/,' MATERIAL PROPERTIES OF INTERBEDS WITH DELAY PROPERTIES'
     1 ,//,'   ZONE        HYDRAULIC           ELASTIC            INEL',
     2 'ASTIC       ',/,'  NUMBER      CONDUCTIVITY     SPECIFIC STORA',
     3 'GE    SPECIFIC STORAGE   ',/,' ',69('-'))
       DO 300 N=1,NMZ
       READ(IN,*) (DP(N,NP),NP=1,3)
  300  CONTINUE
       WRITE(IOUT,305) (N,(DP(N,NP),NP=1,3),N=1,NMZ)
  305  FORMAT(I5,4X,G15.5,5X,G15.5,5X,G15.5)
       LOC3=0
       LOC4=0
       DO 380 KQ=1,NDB
       K=LDN(KQ)
       LOC1=1+(KQ-1)*NCR
C
C18-----READ IN ARRAYS FOR SYSTEMS OF DELAY INTERBEDS.
       IF(IDREST.LE.0) THEN
        WRITE(IOUT,308) KQ
 308    FORMAT(/,1X,' SYSTEM',I4,' OF DELAY BEDS:')
        CALL U2DREL(BUFF(:,:,1),ANAME(5),NROW,NCOL,K,IN,IOUT)
        N1=0
        DO 320 IR=1,NROW
        DO 320 IC=1,NCOL
        N1=N1+1
        LOC2=LOC1+N1-1
        IF(RNB(LOC2).LT.1.0) GO TO 320
        DO 315 N2=1,NN
        LOC3=LOC3+1
        DHP(LOC3)=BUFF(IC,IR,1)
        DH(LOC3)=BUFF(IC,IR,1)
  315   CONTINUE
  320   CONTINUE
        WRITE(IOUT,308) KQ
        CALL U2DREL(BUFF(:,:,1),ANAME(6),NROW,NCOL,K,IN,IOUT)
        N1=0
        DO 330 IR=1,NROW
        DO 330 IC=1,NCOL
        N1=N1+1
        LOC2=LOC1+N1-1
        IF(RNB(LOC2).LT.1.0) GO TO 330
        DO 325 N2=1,NN
        LOC4=LOC4+1
        DHC(LOC4)=BUFF(IC,IR,1)
        IF(DHC(LOC4).GT.DH(LOC4)) DHC(LOC4)=DH(LOC4)
  325   CONTINUE
  330   CONTINUE
       ENDIF
       WRITE(IOUT,308) KQ
C       CALL U2DREL(DCOM(LOC1),ANAME(7),NROW,NCOL,K,IN,IOUT)
       CALL U2DREL(BUFF(:,:,1),ANAME(7),NROW,NCOL,K,IN,IOUT)
       CALL GWF2SUB72D1D(BUFF(:,:,1),NCOL,NROW,DCOM,ND1,LOC1)
       WRITE(IOUT,308) KQ
C       CALL U2DREL(DZ(LOC1),ANAME(8),NROW,NCOL,K,IN,IOUT)
       CALL U2DREL(BUFF(:,:,1),ANAME(8),NROW,NCOL,K,IN,IOUT)
       CALL GWF2SUB72D1D(BUFF(:,:,1),NCOL,NROW,DZ,ND1,LOC1)
       WRITE(IOUT,308) KQ
C       CALL U2DINT(NZ(LOC1),ANAME(9),NROW,NCOL,K,IN,IOUT)
       CALL U2DINT(IBUFF,ANAME(9),NROW,NCOL,K,IN,IOUT)
       L=LOC1-1
       DO 335 I=1,NROW
       DO 335 J=1,NCOL
       L=L+1
       NZ(L)=IBUFF(J,I)
 335   CONTINUE
C
C19-----INITIALIZE ARRAYS FOR SYSTEMS OF DELAY INTERBEDS.
       DO 360 NL=1,NCR
       LOC2=LOC1+NL-1
       IF(RNB(LOC2).GE.1.0.AND.DZ(LOC2).LE.ZERO) THEN
          WRITE(IOUT,355)
 355      FORMAT(' A VALUE OF ZERO WAS FOUND IN THE DZ ARRAY WHERE ',
     1    'DELAY INTERBEDS OCCUR.',/,' MAKE SURE THAT',
     2    ' DZ IS GREATER THAN 0.0 AT ALL CELLS WHERE RNB ',/,
     3    ' IS 1.0 OR MORE. SIMULATION ABORTING')
          CALL USTOP(' ')
       ENDIF
       DZ(LOC2)=DZ(LOC2)*ANNI
  360  CONTINUE
       DO 370 N=1,4
       DVB(KQ,N)=ZERO
  370  CONTINUE
  380  CONTINUE
      ENDIF
C
C20-----SET ALL FLAGS FOR OUTPUT CONTROL TO "FALSE".
      DO 390 I=1,NSTPT
      DO 385 N=1,13
      OCFLGS(N,I)=.FALSE.
  385 CONTINUE
  390 CONTINUE
C
C21-----READ FORMATS AND UNIT NUMBERS OUTPUT FLAGS.
      IF(ISUBOC.GT.0) THEN
       CALL URDCOM(IN,IOUT,LINE)
       LLOC=1
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISBOCF(1),R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISBOCU(1),R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISBOCF(2),R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISBOCU(2),R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISBOCF(3),R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISBOCU(3),R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISBOCF(4),R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISBOCU(4),R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISBOCF(5),R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISBOCU(5),R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISBOCF(6),R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISBOCU(6),R,IOUT,IN)
       WRITE(IOUT,410) (ISBOCF(N),ISBOCU(N),N=1,6)
  410  FORMAT(/,'             SUBSIDENCE PRINT FORMAT IS NUMBER',I4/
     &            '                 UNIT FOR SAVING SUBSIDENCE IS',I4/
     &            '    COMPACTION BY LAYER PRINT FORMAT IS NUMBER',I4/
     &            '        UNIT FOR SAVING COMPACTION BY LAYER IS',I4/
     &            '   COMPACTION BY SYSTEM PRINT FORMAT IS NUMBER',I4/
     &            '       UNIT FOR SAVING COMPACTION BY SYSTEM IS',I4/
     &            '  VERTICAL DISPLACEMENT PRINT FORMAT IS NUMBER',I4/
     &            '      UNIT FOR SAVING VERTICAL DISPLACEMENT IS',I4/
     &            ' NO-DELAY CRITICAL HEAD PRINT FORMAT IS NUMBER',I4/
     &            '     UNIT FOR SAVING NO-DELAY CRITICAL HEAD IS',I4/
     &            '    DELAY CRITICAL HEAD PRINT FORMAT IS NUMBER',I4/
     &            '        UNIT FOR SAVING DELAY CRITICAL HEAD IS',I4)
       NTSSUM(1)=0
       IF(NPER.GT.1) THEN
        DO 415 N=2,NPER
        NTSSUM(N)=NTSSUM(N-1)+NSTP(N-1)
  415   CONTINUE
       ENDIF
       DO 450 NOCLIN=1,ISUBOC
       CALL URDCOM(IN,IOUT,LINE)
       LLOC=1
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISP1,R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,ISP2,R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,JTS1,R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,JTS2,R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IFL(1),R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IFL(2),R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IFL(3),R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IFL(4),R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IFL(5),R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IFL(6),R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IFL(7),R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IFL(8),R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IFL(9),R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IFL(10),R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IFL(11),R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IFL(12),R,IOUT,IN)
       CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IFL(13),R,IOUT,IN)
       IF(ISP1.LT.1) ISP1=1
       IF(ISP1.GT.NPER) ISP1=NPER
       IF(ISP2.LT.1) ISP2=1
       IF(ISP2.GT.NPER) ISP2=NPER
       IF(ISP1.GT.ISP2) ISP1=ISP2
       DO 440 I=ISP1,ISP2
       J1=JTS1
       J2=JTS2
       IF(J1.LT.1) J1=1
       IF(J1.GT.NSTP(I)) J1=NSTP(I)
       IF(J2.LT.1) J2=1
       IF(J2.GT.NSTP(I)) J2=NSTP(I)
       IF(J1.GT.J2) J1=J2
       DO 430 J=J1,J2
       ILOC=NTSSUM(I)+J
       DO 420 N=1,13
       IF(IFL(N).GT.0) OCFLGS(N,ILOC)=.TRUE.
       IF(IFL(N).EQ.0) OCFLGS(N,ILOC)=.FALSE.
  420  CONTINUE
  430  CONTINUE
  440  CONTINUE
  450  CONTINUE
      ENDIF
C
C22-----RETURN
  500 CALL SGWF2SUB7PSV(IGRID)
      RETURN
      END
      SUBROUTINE GWF2SUB7ST(KPER,IGRID)
C     ******************************************************************
C        SET PRECONSOLIDATION HEAD (HC AND DHC) EQUAL TO THE STEADY-
C        STATE HEAD IF HEAD IS LOWER THAN PRECONSOLIDATION HEAD.
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,    ONLY:HNEW,NCOL,NROW,ISSFLG
      USE GWFSUBMODULE ,ONLY: RNB,LN,LDN,HC,DHP,DH,DHC,NZ,DZ,
     1                        NDF,NNDF,NNDB,NDB,NN
C
C     ------------------------------------------------------------------
      CALL SGWF2SUB7PNT(IGRID)
C
C1------RETURN IF THIS IS NOT THE SECOND STRESS PERIOD OR IF THE FIRST
C1------STRESS PERIOD WAS TRANSIENT.
      IF(KPER.NE.2) RETURN
      IF(ISSFLG(1).EQ.0) RETURN
C2-----MAKE SURE THAT NO-DELAY PRECONSOLIDATION HEAD VALUES ARE CONSISTENT
C2-----WITH STEADY-STATE HEADS.
      NCR=NROW*NCOL
      IF(NNDF) THEN
       DO 20 KQ=1,NNDB
       K=LN(KQ)
       NQ=(KQ-1)*NCR
C       NK=(K-1)*NCR
       DO 10 IR=1,NROW
       NQR=NQ+(IR-1)*NCOL
C       NKR=NK+(IR-1)*NCOL
       DO 10 IC=1,NCOL
       LOC2=NQR+IC
C       LOC2H=NKR+IC
       HTMP=HNEW(IC,IR,K)
       IF(HC(LOC2).GT.HTMP) HC(LOC2)=HTMP
   10  CONTINUE
   20  CONTINUE
      ENDIF
C3-----MAKE SURE THAT DELAY PRECONSOLIDATION HEAD VALUES ARE CONSISTENT
C3-----WITH STEADY-STATE HEADS. ALSO, SET HEAD (DH) AND HEAD FOR
C3-----PREVIOUS TIME STEP (DHP) EQUAL TO THE AQUIFER HEAD (HNEW).
      IF(NDF) THEN
       LOC3=0
       DO 40 KQ=1,NDB
       K=LDN(KQ)
       NQ=(KQ-1)*NCR
C       NK=(K-1)*NCR
       N1=0
       DO 30 IR=1,NROW
       DO 30 IC=1,NCOL
       N1=N1+1
       LOC2=NQ+N1
C       LOC2H=NK+N1
       HTMP=HNEW(IC,IR,K)
       IF(RNB(LOC2).LT.1.0) GO TO 30
       DO 18 N2=1,NN
       LOC3=LOC3+1
       IF(DHC(LOC3).GT.HTMP) DHC(LOC3)=HTMP
       DH(LOC3)=HTMP
       DHP(LOC3)=HTMP
   18  CONTINUE
   30  CONTINUE
   40  CONTINUE
      ENDIF
C4-----RETURN.
      RETURN
      END
      SUBROUTINE GWF2SUB7FM(KPER,KITER,ISIP,IGRID)
C     ******************************************************************
C        ADD INTERBED STORAGE TERMS TO RHS AND HCOF
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,    ONLY: RHS,HCOF,HNEW,HOLD,IBOUND,DELR,DELC,
     1                     NCOL,NROW,ISSFLG
      USE GWFBASMODULE, ONLY: DELT
      USE SIPMODULE,    ONLY: V,HCLOSE
      USE GWFSUBMODULE ,ONLY: RNB,LN,LDN,HC,SCE,SCV,DHP,DH,DHC,NZ,DZ,DP,
     1                        BB,NDF,NNDF,AC1,AC2,ITMIN,NN,
     1                        NDB,NNDB,NMZ
      LOGICAL ICHK
C
      DATA ICHK/.FALSE./
C     ------------------------------------------------------------------
      CALL SGWF2SUB7PNT(IGRID)
      ZERO=0.0
C
C0------SKIP CALCULATIONS IF THIS IS A STEADY-STATE STRESS PERIOD.
      IF(ISSFLG(KPER).EQ.1) RETURN
C
C1------INITIALIZE
      TLED=1./DELT
      NCR=NCOL*NROW
C
      IF(NNDF) THEN
C2------FIND LAYERS WITH INTERBED STORAGE
       DO 110 KQ=1,NNDB
       K=LN(KQ)
       LOCT=(KQ-1)*NCR
       N=0
       DO 100 IR=1,NROW
       DO 100 IC=1,NCOL
       N=N+1
       LOC2=LOCT+N
       IF(IBOUND(IC,IR,K).LE.0) GO TO 100
C
C3------DETERMINE STORAGE CAPACITIES FOR CELL AT START AND END OF STEP
       RHO1=SCE(LOC2)*TLED
       RHO2=RHO1
       HCTMP=HC(LOC2)
       IF(HNEW(IC,IR,K).LT.HCTMP) RHO2=SCV(LOC2)*TLED
C
C4------ADD APPROPRIATE TERMS TO RHS AND HCOF
       RHS(IC,IR,K)=RHS(IC,IR,K)-HCTMP*(RHO2-RHO1)-RHO1*HOLD(IC,IR,K)
       HCOF(IC,IR,K)=HCOF(IC,IR,K)-RHO2
  100  CONTINUE
  110  CONTINUE
      ENDIF
      IF(NDF) THEN
       LOC3=1-NN
       DO 420 KQ=1,NDB
       K=LDN(KQ)
       NQ=(KQ-1)*NCR
       DO 410 IR=1,NROW
       NQR=NQ+(IR-1)*NCOL
       DO 410 IC=1,NCOL
       LOC2=NQR+IC
       RNB2=RNB(LOC2)
       IF(RNB2.LT.1.0) GO TO 410
       LOC3=LOC3+NN
       IF(IBOUND(IC,IR,K).LE.0) GO TO 410
       IF(ISIP.GT.0.AND.ICHK) THEN
        VV=V(IC,IR,K)
       ELSE
        VV=ZERO
        ICHK=.TRUE.
       ENDIF
       NZONE=NZ(LOC2)
       DZZ=DZ(LOC2)
       CI=DP(NZONE,1)/DZZ
       IF(ISIP.GT.0) THEN
         IF(KITER.GT.ITMIN.AND.ABS(VV).LT.HCLOSE) GO TO 205
       END IF
C
C5------ASSEMBLE COEFFICIENTS FOR DIRECT SOLUTION OF HEAD IN INTERBED
       HAQ=HNEW(IC,IR,K)+VV*AC1
       SSE=DP(NZONE,2)
       SSV=DP(NZONE,3)
       NEND=LOC3+NN-1
       CALL SGWF2SUB7A(HAQ,TLED,CI,SSE,SSV,DZZ,
     1      DH(LOC3:NEND),DHP(LOC3:NEND),DHC(LOC3:NEND),NN)
C
C6------SOLVE FOR HEAD CHANGES IN STRING USING GAUSSIAN ELIMINATION.
C6------ADD CHANGES TO HEAD VALUES TO GET HEAD AT CURRENT ITERATION.
       CALL SGWF2SUB7S(NN)
       DO 200 N=1,NN
       DH(LOC3+N-1)=DH(LOC3+N-1)+BB(N)*AC2
  200  CONTINUE
C
C7------CALCULATE STORAGE CHANGE IN INTERBEDS
  205  AREA=DELR(IC)*DELC(IR)
       STRGS=ZERO
       L1=LOC3
       L2=LOC3+NN-1
       DO 210 LOC4=L1,L2
       HHOLD=DHP(LOC4)
       HHNEW=DH(LOC4)
       HHC=DHC(LOC4)
C
C8------GET STORAGE CAPACITIES AT BEGINNING AND END OF TIME STEP.
       SBGN=DP(NZ(LOC2),2)
       SEND=SBGN
       IF(HHNEW.LT.HHC) SEND=DP(NZ(LOC2),3)
C
C9------CALCULATE VOLUME CHANGE IN INTERBED STORAGE FOR TIME STEP.
       STR1=(HHC*(SEND-SBGN)+SBGN*HHOLD-
     1                 SEND*HHNEW)*DZ(LOC2)*RNB(LOC2)*2.
       IF(LOC4.EQ.L2) STR1=STR1*.5
       STRGS=STRGS+STR1
  210  CONTINUE
       RATES=STRGS*AREA*TLED
C
C10-----ADD APPROPRIATE TERMS TO RHS AND HCOF
       RHS(IC,IR,K)=RHS(IC,IR,K)-RATES
  410  CONTINUE
  420  CONTINUE
      ENDIF
C
C11-----RETURN
      RETURN
      END
      SUBROUTINE GWF2SUB7BD(KSTP,KPER,IGRID)
C     ******************************************************************
C     CALCULATE VOLUMETRIC BUDGET FOR INTERBED STORAGE
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,       ONLY: IOUT,NCOL,NROW,NLAY,IBOUND,HNEW,HOLD,
     1                        BUFF,DELR,DELC,ISSFLG
      USE GWFBASMODULE, ONLY: VBVL,VBNM,MSUM,ICBCFL,DELT
      USE GWFSUBMODULE ,ONLY: RNB,LN,LDN,HC,SCE,SCV,SUB,DHP,DH,DHC,
     1                        NZ,DZ,DCOM,DP,DVB,NDF,NNDF,NN,ND2,NDB,
     2                        NNDB,NMZ,IIBSCB
      CHARACTER*16 TEXT(2)
C
      DATA TEXT(1) /'INST. IB STORAGE'/
      DATA TEXT(2) /'DELAY IB STORAGE'/
C     ------------------------------------------------------------------
      CALL SGWF2SUB7PNT(IGRID)
      ZERO=0.0
      TLED=ZERO
      IF(ISSFLG(KPER).EQ.0) TLED=1./DELT
C
C1------INITIALIZE CELL-BY-CELL FLOW TERM FLAG (IBD) AND
C1------SET IF CELL-BY-CELL FLOW TERMS ARE NEEDED.
      IBD=0
      IF(ICBCFL.NE.0 .AND. IIBSCB.GT.0 ) IBD=1
      NCR=NCOL*NROW
C
C2------RUN THROUGH EVERY CELL IN THE GRID WITH INTERBED STORAGE.
      IF(NNDF) THEN
C
C3-------CELL-BY-CELL FLOW TERMS ARE NEEDED SET IBD AND CLEAR BUFFER.
       IF(IBD.EQ.1) THEN
        DO 90 K=1,NLAY
        DO 90 IR=1,NROW
        DO 90 IC=1,NCOL
        BUFF(IC,IR,K)=ZERO
   90   CONTINUE
       ENDIF
       STOIN=ZERO
       STOUT=ZERO
C
C4------IF THIS IS A STEADY-STATE STRESS PERIOD, SKIP CALCULATIONS
       IF(ISSFLG(KPER).EQ.1) GO TO 111
C
C5------CALCULATE NO-DELAY INTERBED STORAGE CHANGE FOR EACH CELL
       DO 110 KQ=1,NNDB
       K=LN(KQ)
       NQ=(KQ-1)*NCR
       DO 100 IR=1,NROW
       NQR=NQ+(IR-1)*NCOL
       DO 100 IC=1,NCOL
       LOC2=NQR+IC
C
C6------CALCULATE FLOW FROM STORAGE (VARIABLE HEAD CELLS ONLY)
       IF(IBOUND(IC,IR,K).LE.0) GO TO 100
       HHOLD=HOLD(IC,IR,K)
       HHNEW=HNEW(IC,IR,K)
       HHC=HC(LOC2)
C
C7------GET STORAGE CAPACITIES AT BEGINNING AND END OF TIME STEP.
       SBGN=SCE(LOC2)
       SEND=SBGN
       IF(HHNEW.LT.HHC) SEND=SCV(LOC2)
C
C8------CALCULATE VOLUME CHANGE IN INTERBED STORAGE FOR TIME STEP.
       STRG=HHC*(SEND-SBGN)+SBGN*HHOLD-SEND*HHNEW
C
C9------ACCUMULATE SUBSIDENCE ASSOCIATED WITH CHANGE IN STORAGE
       SUB(LOC2)=SUB(LOC2)+STRG/(DELR(IC)*DELC(IR))
C
C10-----IF C-B-C FLOW TERMS ARE TO BE SAVED THEN ADD RATE TO BUFFER.
       IF(IBD.EQ.1) BUFF(IC,IR,K)=BUFF(IC,IR,K)+STRG*TLED
C
C11-----SEE IF FLOW IS INTO OR OUT OF STORAGE.
       IF(STRG.LE.ZERO) THEN
        STOUT=STOUT-STRG
       ELSE
        STOIN=STOIN+STRG
       ENDIF
  100  CONTINUE
  110  CONTINUE
C
C12-----IF C-B-C FLOW TERMS WILL BE SAVED CALL UBUDSV TO RECORD THEM.
  111  IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT(1),IIBSCB,BUFF,NCOL,
     1                           NROW,NLAY,IOUT)
C
C13-----MOVE RATES,VOLUMES & LABELS INTO ARRAYS FOR PRINTING.
       VBVL(3,MSUM)=STOIN*TLED
       VBVL(4,MSUM)=STOUT*TLED
       VBVL(1,MSUM)=VBVL(1,MSUM)+STOIN
       VBVL(2,MSUM)=VBVL(2,MSUM)+STOUT
       VBNM(MSUM)=TEXT(1)
C
C14-----INCREMENT BUDGET TERM COUNTER
       MSUM=MSUM+1
C
C15-----UPDATE PRECONSOLIDATION HEAD ARRAY
       DO 310 KQ=1,NNDB
       K=LN(KQ)
       LOCT=(KQ-1)*NCR
       N=0
       DO 300 IR=1,NROW
       DO 300 IC=1,NCOL
       N=N+1
       LOC2=LOCT+N
       IF(IBOUND(IC,IR,K).LE.0) GO TO 300
       HHNEW=HNEW(IC,IR,K)
       IF(HHNEW.LT.HC(LOC2)) HC(LOC2)=HHNEW
  300  CONTINUE
  310  CONTINUE
      ENDIF
      IF(NDF) THEN
       IF(IBD.EQ.1) THEN
        DO 390 K=1,NLAY
        DO 390 IR=1,NROW
        DO 390 IC=1,NCOL
        BUFF(IC,IR,K)=ZERO
  390   CONTINUE
       ENDIF
       RATIN=ZERO
       RATOUT=ZERO
C
C16-----IF THIS IS A STEADY-STATE STRESS PERIOD, SKIP CALCULATIONS
       IF(ISSFLG(KPER).EQ.1) GO TO 421
C
C17-----CALCULATE NO-DELAY INTERBED STORAGE CHANGE FOR EACH CELL
       LOC3=1-NN
       DO 420 KQ=1,NDB
       K=LDN(KQ)
       NQ=(KQ-1)*NCR
       STRGT=ZERO
       RATBSM=ZERO
       DO 410 IR=1,NROW
       NQR=NQ+(IR-1)*NCOL
       DO 410 IC=1,NCOL
       AREA=DELR(IC)*DELC(IR)
       LOC2=NQR+IC
       RNB2=RNB(LOC2)
       IF(RNB2.LT.1.0) GO TO 410
       LOC3=LOC3+NN
       IF(IBOUND(IC,IR,K).LE.0) GO TO 410
       HD1=DH(LOC3)
       HHNEW=HNEW(IC,IR,K)
C
C18-----CALCULATE CONDUCTANCE BETWEEN AQUIFER AND FIRST CELL IN INTERBED
C18-----ACCOUNTING FOR BOTH HALVES OF RNB(LOC2) BEDS IN SYSTEM
       COND=4.*RNB2*DP(NZ(LOC2),1)*AREA/DZ(LOC2)
C
C19-----CALCULATE THE FLOW RATE INTO THE CELL
       RATB=DBLE(COND)*(DBLE(HD1)-HNEW(IC,IR,K))
       STRGS=ZERO
       L1=LOC3
       L2=LOC3+NN-1
       DO 401 LOC4=L1,L2
       HHOLD=DHP(LOC4)
       HHNEW=DH(LOC4)
       HHC=DHC(LOC4)
C
C20-----GET STORAGE CAPACITIES AT BEGINNING AND END OF TIME STEP.
       SBGN=DP(NZ(LOC2),2)
       SEND=SBGN
       IF(HHNEW.LT.HHC) SEND=DP(NZ(LOC2),3)
C
C21-----CALCULATE VOLUME CHANGE IN INTERBED STORAGE FOR TIME STEP.
       STR1=(HHC*(SEND-SBGN)+SBGN*HHOLD-
     1                        SEND*HHNEW)*DZ(LOC2)*RNB(LOC2)*2.
       IF(LOC4.EQ.L2) STR1=STR1*.5
       STRGS=STRGS+STR1
C
C22-----ACCUMULATE SUBSIDENCE ASSOCIATED WITH CHANGE IN STORAGE
       DCOM(LOC2)=DCOM(LOC2)+STR1
  401  CONTINUE
       STRGS=STRGS*AREA
       STRGT=STRGT+STRGS
       RATS=STRGS*TLED
       IF(IBD.EQ.1) BUFF(IC,IR,K)=BUFF(IC,IR,K)+RATS
       RATBSM=RATBSM-RATS
       IF(RATS.LE.ZERO) THEN
        RATOUT=RATOUT-RATS
       ELSE
        RATIN=RATIN+RATS
       ENDIF
  410  CONTINUE
       DVB(KQ,1)=DVB(KQ,1)+STRGT
       DVB(KQ,2)=DVB(KQ,2)+RATBSM*DELT
       DVB(KQ,3)=STRGT*TLED
       DVB(KQ,4)=RATBSM
  420  CONTINUE
C
C23-----IF C-B-C FLOW TERMS WILL BE SAVED CALL UBUDSV TO RECORD THEM.
  421  IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT(2),IIBSCB,BUFF,
     1                           NCOL,NROW,NLAY,IOUT)
C
C24-----MOVE RATES,VOLUMES & LABELS INTO ARRAYS FOR PRINTING.
       VBVL(3,MSUM)=RATIN
       VBVL(4,MSUM)=RATOUT
       VBVL(1,MSUM)=VBVL(1,MSUM)+RATIN*DELT
       VBVL(2,MSUM)=VBVL(2,MSUM)+RATOUT*DELT
       VBNM(MSUM)=TEXT(2)
C
C25-----INCREMENT BUDGET TERM COUNTER
       MSUM=MSUM+1
C
C26-----UPDATE PRECONSOLIDATION HEAD ARRAY, PREVIOUS HEAD ARRAY FOR
C26-----SYSTEMS OF DELAY INTERBEDS.
       DO 500 N=1,ND2
       IF(DH(N).LT.DHC(N)) DHC(N)=DH(N)
       DHP(N)=DH(N)
  500  CONTINUE
      ENDIF
C
C27----RETURN
      RETURN
      END
      SUBROUTINE GWF2SUB7OT(KSTP,KPER,IN,IGRID)
C     ******************************************************************
C     PRINT AND STORE SUBSIDENCE, COMPACTION AND CRITICAL HEAD.
C     ******************************************************************
C
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,     ONLY: IOUT,NCOL,NROW,NLAY,NSTP,BUFF,ISSFLG
      USE GWFBASMODULE, ONLY:PERTIM,TOTIM
      USE GWFSUBMODULE ,ONLY: LN,LDN,SUB,HC,RNB,DCOM,DHC,DVB,NDF,NNDF,
     &                        NTSSUM,OCFLGS,OCLAY,ILSYS,ISBOCF,ISBOCU,
     &                        NN,NNDB,NDB
      CHARACTER*16 TEXT(7)
      LOGICAL IBDPR
      DATA TEXT
     &  /'      SUBSIDENCE',      'LAYER COMPACTION',
     &   'NDSYS COMPACTION',      ' DSYS COMPACTION',
     &   '  Z DISPLACEMENT',      'ND CRITICAL HEAD',
     2   ' D CRITICAL HEAD'/
C     ------------------------------------------------------------------
      CALL SGWF2SUB7PNT(IGRID)
C
      ZERO=0.0
      NCR=NCOL*NROW
C1------INITIALIZE TIME STEP POINTER TO RETRIEVE FLAGS FOR PRINTING AND
C1------SAVING ARRAYS. SET FLAG FOR PRINTING BUDGET FOR DELAY INTERBEDS.
      NNSTP=NTSSUM(KPER)+KSTP
      IBDPR=.FALSE.
      IF(KSTP.EQ.NSTP(KPER).OR.OCFLGS(13,NNSTP)) IBDPR=.TRUE.
      IF(ISSFLG(KPER).EQ.1) IBDPR=.FALSE.
C
C3------PRINT AND STORE SUBSIDENCE, FIRST, CLEAR OUT BUFF.
      IF(OCFLGS(1,NNSTP).OR.OCFLGS(2,NNSTP)) THEN
       DO 30 K=1,NLAY
       DO 30 IR=1,NROW
       DO 30 IC=1,NCOL
       BUFF(IC,IR,K)=ZERO
   30  CONTINUE
C
C4-------SUM COMPACTION IN ALL LAYERS TO GET SUBSIDENCE.
       IF(NNDF) THEN
        DO 50 KQ=1,NNDB
        LOCT=(KQ-1)*NCR
        N=0
        DO 40 IR=1,NROW
        DO 40 IC=1,NCOL
        N=N+1
        LOC2=LOCT+N
        BUFF(IC,IR,1)=BUFF(IC,IR,1)+SUB(LOC2)
   40   CONTINUE
   50   CONTINUE
       ENDIF
       IF(NDF) THEN
        DO 70 KQ=1,NDB
        LOCT=(KQ-1)*NCR
        N=0
        DO 60 IR=1,NROW
        DO 60 IC=1,NCOL
        N=N+1
        LOC2=LOCT+N
        BUFF(IC,IR,1)=BUFF(IC,IR,1)+DCOM(LOC2)
   60   CONTINUE
   70   CONTINUE
       ENDIF
C
C5-------PRINT SUBSIDENCE.
       IF(OCFLGS(1,NNSTP)) THEN
        IF(ISBOCF(1).LT.0) CALL ULAPRS(BUFF(:,:,1),TEXT(1),KSTP,KPER,
     1            NCOL,NROW,1,-ISBOCF(1),IOUT)
        IF(ISBOCF(1).GE.0) CALL ULAPRW(BUFF(:,:,1),TEXT(1),KSTP,KPER,
     1            NCOL,NROW,1,ISBOCF(1),IOUT)
       ENDIF
C
C6-------STORE SUBSIDENCE.
       IF(OCFLGS(2,NNSTP)) THEN
        CALL ULASAV(BUFF(:,:,1),TEXT(1),KSTP,KPER,PERTIM,TOTIM,NCOL,
     1              NROW,1,ISBOCU(1))
       ENDIF
      ENDIF
C
C7------PRINT AND STORE COMPACTION FOR EACH SYSTEM OF INTERBEDS,
C7------INCLUDING DELAY AND NO-DELAY INTERBEDS.
      IF(OCFLGS(5,NNSTP).OR.OCFLGS(6,NNSTP)) THEN
       IF(NNDF) THEN
        DO 80 KQ=1,NNDB
        K=LN(KQ)
        LOC2=(KQ-1)*NCR+1
        IF(OCFLGS(5,NNSTP)) THEN
       WRITE(IOUT,76) KQ
  76   FORMAT(/,1X,' SYSTEM',I4,' OF NO-DELAY BEDS:')
        NEND=LOC2+NCOL*NROW-1
        IF(ISBOCF(3).LT.0) CALL ULAPRS(SUB(LOC2:NEND),TEXT(3),KSTP,
     1            KPER,NCOL,NROW,K,-ISBOCF(3),IOUT)
        IF(ISBOCF(3).GE.0) CALL ULAPRW(SUB(LOC2:NEND),TEXT(3),KSTP,
     1            KPER,NCOL,NROW,K,ISBOCF(3),IOUT)
        ENDIF
        IF(OCFLGS(6,NNSTP)) THEN
         CALL ULASAV(SUB(LOC2:NEND),TEXT(3),KSTP,KPER,PERTIM,TOTIM,
     1              NCOL,NROW,KQ,ISBOCU(3))
        ENDIF
   80   CONTINUE
       ENDIF
       IF(NDF) THEN
        DO 90 KQ=1,NDB
        K=LDN(KQ)
        LOC2=(KQ-1)*NCR+1
        IF(OCFLGS(5,NNSTP)) THEN
        WRITE(IOUT,82) KQ
  82   FORMAT(/,1X,' SYSTEM',I4,' OF DELAY BEDS:')
         NEND=LOC2+NCOL*NROW-1
         IF(ISBOCF(3).LT.0) CALL ULAPRS(DCOM(LOC2:NEND),TEXT(4),KSTP,
     1            KPER,NCOL,NROW,K,-ISBOCF(3),IOUT)
         IF(ISBOCF(3).GE.0) CALL ULAPRW(DCOM(LOC2:NEND),TEXT(4),KSTP,
     1            KPER,NCOL,NROW,K,ISBOCF(3),IOUT)
        ENDIF
        IF(OCFLGS(6,NNSTP)) THEN
         CALL ULASAV(DCOM(LOC2:NEND),TEXT(4),KSTP,KPER,PERTIM,TOTIM,
     1              NCOL,NROW,KQ,ISBOCU(3))
        ENDIF
   90   CONTINUE
       ENDIF
      ENDIF
C
C8------SUM COMPACTION IN EACH LAYER IN THE BUFF ARRAY FOR SAVING
C8------OR PRINTING COMPACTION OR VERTICAL DISPLACEMENT BY MODEL
C8------LAYER. FIRST, CLEAR OUT BUFF.
      IF(OCFLGS(3,NNSTP).OR.OCFLGS(4,NNSTP).OR.
     & OCFLGS(7,NNSTP).OR.OCFLGS(8,NNSTP)) THEN
       DO 93 NL=1,NLAY
       OCLAY(NL)=.FALSE.
   93  CONTINUE
       DO 95 K=1,NLAY
       DO 95 IR=1,NROW
       DO 95 IC=1,NCOL
       BUFF(IC,IR,K)=ZERO
   95  CONTINUE
C
C9-------SUM NO-DELAY COMPACTION IN ALL MODEL LAYERS.
       IF(NNDF) THEN
        DO 105 KQ=1,NNDB
        K=LN(KQ)
        OCLAY(K)=.TRUE.
        LOCT=(KQ-1)*NCR
        N=0
        DO 100 IR=1,NROW
        DO 100 IC=1,NCOL
        N=N+1
        LOC2=LOCT+N
        BUFF(IC,IR,K)=BUFF(IC,IR,K)+SUB(LOC2)
  100   CONTINUE
  105   CONTINUE
       ENDIF
C
C10------SUM DELAY COMPACTION IN ALL MODEL LAYERS.
       IF(NDF) THEN
        DO 115 KQ=1,NDB
        K=LDN(KQ)
        OCLAY(K)=.TRUE.
        LOCT=(KQ-1)*NCR
        N=0
        DO 110 IR=1,NROW
        DO 110 IC=1,NCOL
        N=N+1
        LOC2=LOCT+N
        BUFF(IC,IR,K)=BUFF(IC,IR,K)+DCOM(LOC2)
  110   CONTINUE
  115   CONTINUE
       ENDIF
      ENDIF
C
C11-----PRINT COMPACTION BY LAYER.
      IF(OCFLGS(3,NNSTP)) THEN
       DO 120 KL=1,NLAY
       IF(.NOT.OCLAY(KL)) GO TO 120
       KKL=KL
       IF(ISBOCF(2).LT.0) CALL ULAPRS(BUFF(:,:,KL),TEXT(2),KSTP,KPER,
     1           NCOL,NROW,KKL,-ISBOCF(2),IOUT)
       IF(ISBOCF(2).GE.0) CALL ULAPRW(BUFF(:,:,KL),TEXT(2),KSTP,KPER,
     1            NCOL,NROW,KKL,ISBOCF(2),IOUT)
  120  CONTINUE
      ENDIF
C
C12-----STORE COMPACTION BY LAYER.
      IF(OCFLGS(4,NNSTP)) THEN
       DO 125 KL=1,NLAY
       IF(.NOT.OCLAY(KL)) GO TO 125
       KKL=KL
       CALL ULASAV(BUFF(:,:,KL),TEXT(2),KSTP,KPER,PERTIM,TOTIM,NCOL,
     1             NROW,KKL,ISBOCU(2))
  125  CONTINUE
      ENDIF
C
C13----CALCULATE VERTICAL DISPLACEMENT.
      IF(OCFLGS(7,NNSTP).OR.OCFLGS(8,NNSTP)) THEN
       NL1=NLAY-1
       IF(NLAY.GT.1) THEN
        DO 140 KL=NL1,1,-1
        LOCT3=(KL-1)*NCR
        N=0
        DO 135 IR=1,NROW
        DO 135 IC=1,NCOL
        N=N+1
        BUFF(IC,IR,KL)=BUFF(IC,IR,KL)+BUFF(IC,IR,KL+1)
  135   CONTINUE
  140   CONTINUE
       ENDIF
C
C14-----PRINT VERTICAL DISPLACEMENT FOR ALL MODEL LAYERS.
       IF(OCFLGS(7,NNSTP)) THEN
        DO 145 KL=1,NLAY
        KKL=KL
        IF(ISBOCF(4).LT.0) CALL ULAPRS(BUFF(:,:,KL),TEXT(5),KSTP,KPER,
     1            NCOL,NROW,KKL,-ISBOCF(4),IOUT)
        IF(ISBOCF(4).GE.0) CALL ULAPRW(BUFF(:,:,KL),TEXT(5),KSTP,KPER,
     1             NCOL,NROW,KKL,ISBOCF(4),IOUT)
  145   CONTINUE
       ENDIF
C
C15-----SAVE VERTICAL DISPLACEMENT FOR ALL MODEL LAYERS.
       IF(OCFLGS(8,NNSTP)) THEN
        DO 150 KL=1,NLAY
        KKL=KL
        CALL ULASAV(BUFF(:,:,KL),TEXT(5),KSTP,KPER,PERTIM,TOTIM,NCOL,
     1              NROW,KKL,ISBOCU(4))
  150   CONTINUE
       ENDIF
      ENDIF
C
C16-----PRINT CRITICAL HEAD FOR SYSTEMS OF NO-DELAY INTERBEDS.
C16-----STORAGE.
      IF(OCFLGS(9,NNSTP).OR.OCFLGS(10,NNSTP)) THEN
       IF(NNDF) THEN
       DO 155 NL=1,NLAY
       OCLAY(NL)=.TRUE.
  155  CONTINUE
        DO 160 KQ=1,NNDB
        K=LN(KQ)
        LOC2=(KQ-1)*NCR+1
        IF(.NOT.OCLAY(K)) GO TO 160
        NSYS=1
        ILSYS(1)=KQ
        IF(KQ.LT.NNDB) THEN
         KS1=KQ+1
         DO 152 KS=KS1,NNDB
         IF(LN(KS).EQ.K) THEN
          NSYS=NSYS+1
          ILSYS(NSYS)=KS
         ENDIF
  152    CONTINUE
        ENDIF
        OCLAY(K)=.FALSE.
        IF(OCFLGS(9,NNSTP)) THEN
         WRITE(IOUT,154) (ILSYS(NS),NS=1,NSYS)
  154  FORMAT(/,1X,' SYSTEM OR SYSTEMS OF NO-DELAY BEDS:',20I3)
         NEND=LOC2+NCOL*NROW-1
         IF(ISBOCF(5).LT.0) CALL ULAPRS(HC(LOC2:NEND),TEXT(6),KSTP,KPER,
     1            NCOL,NROW,K,-ISBOCF(5),IOUT)
         IF(ISBOCF(5).GE.0) CALL ULAPRW(HC(LOC2:NEND),TEXT(6),KSTP,KPER,
     1             NCOL,NROW,K,ISBOCF(5),IOUT)
        ENDIF
        IF(OCFLGS(10,NNSTP)) THEN
         CALL ULASAV(HC(LOC2:NEND),TEXT(6),KSTP,KPER,PERTIM,TOTIM,NCOL,
     1            NROW,K,ISBOCU(5))
        ENDIF
  160   CONTINUE
       ENDIF
      ENDIF
C
C17-----PRINT CRITICAL HEAD FOR ALL SYSTEMS OF DELAY INTERBED.
      IF(OCFLGS(11,NNSTP).OR.OCFLGS(12,NNSTP)) THEN
       IF(NDF) THEN
        LOC4=0
        DO 190 KQ=1,NDB
        K=LDN(KQ)
        LOCT=(KQ-1)*NCR
        N=0
        DO 180 IR=1,NROW
        DO 180 IC=1,NCOL
        N=N+1
        BUFF(IC,IR,1)=ZERO
        LOC2=LOCT+N
        IF(RNB(LOC2).LT.1.0) GO TO 180
        LOC4=LOC4+NN
        BUFF(IC,IR,1)=DHC(LOC4)
  180   CONTINUE
        IF(OCFLGS(11,NNSTP)) THEN
         WRITE(IOUT,82) KQ
         IF(ISBOCF(6).LT.0) CALL ULAPRS(BUFF(:,:,1),TEXT(7),KSTP,KPER,
     1             NCOL,NROW,K,-ISBOCF(6),IOUT)
         IF(ISBOCF(6).GE.0) CALL ULAPRW(BUFF(:,:,1),TEXT(7),KSTP,KPER,
     1              NCOL,NROW,K,ISBOCF(6),IOUT)
        ENDIF
        IF(OCFLGS(12,NNSTP)) THEN
         CALL ULASAV(BUFF(:,:,1),TEXT(7),KSTP,KPER,PERTIM,TOTIM,NCOL,
     1            NROW,KQ,ISBOCU(6))
        ENDIF
  190   CONTINUE
       ENDIF
      ENDIF
C
C18-----PRINT VOLUMETRIC BUDGET FOR SYSTEMS OF DELAY INTERBEDS
      IF(NDF.AND.IBDPR) THEN
       WRITE(IOUT,230) KSTP,KPER
       SUMSV=ZERO
       SUMBV=ZERO
       SUMSR=ZERO
       SUMBR=ZERO
       DO 200 KQ=1,NDB
       DISCV=ZERO
       DISCR=ZERO
       SUMV=DVB(KQ,1)+DVB(KQ,2)
       SUMR=DVB(KQ,3)+DVB(KQ,4)
       SUMSV=SUMSV+DVB(KQ,1)
       SUMBV=SUMBV+DVB(KQ,2)
       SUMSR=SUMSR+DVB(KQ,3)
       SUMBR=SUMBR+DVB(KQ,4)
       IF(DVB(KQ,1).NE.ZERO) DISCV=100.*SUMV/DVB(KQ,1)
       IF(DVB(KQ,3).NE.ZERO) DISCR=100.*SUMR/DVB(KQ,3)
       WRITE(IOUT,240) KQ,DVB(KQ,1),DVB(KQ,2),SUMV,DISCV,
     1                    DVB(KQ,3),DVB(KQ,4),SUMR,DISCR
  200  CONTINUE
       IF(NDB.GT.1) THEN
        DISCV=ZERO
        DISCR=ZERO
        SUMV=SUMSV+SUMBV
        SUMR=SUMSR+SUMBR
        IF(SUMSV.NE.ZERO) DISCV=100.*SUMV/SUMSV
        IF(SUMSR.NE.ZERO) DISCR=100.*SUMR/SUMSR
        WRITE(IOUT,250) SUMSV,SUMBV,SUMV,DISCV,
     1                  SUMSR,SUMBR,SUMR,DISCR
       ENDIF
      ENDIF
C
C19-----RETURN
      RETURN
C
C20-----FORMATS
C
  230 FORMAT(/,31X,'VOLUMETRIC BUDGET FOR SYSTEMS OF DELAY INTERBEDS',
     1 /,42X,'AT END OF TIME STEP',I3,' IN ','STRESS PERIOD',I3,
     2 //,'          |   C U M U L A T I V E   ',
     3 'V O L U M E S   L**3           | R A T E S   F O R  T H I S',
     4 '  T I M E  S T E P   L**3/T  |',/,'   SYSTEM |    CHANGE IN',
     5 '     BOUNDARY                    PERCENT   |    CHANGE IN  ',
     6 '   BOUNDARY                    PERCENT   |',/,'   NUMBER | ',
     7 '    STORAGE        FLOW           SUM      DISCREPANCY |   ',
     8 '  STORAGE        FLOW           SUM      DISCREPANCY |',/,2X,
     9 8('-'),'|',56('-'),'|',56('-'),'|')
  240 FORMAT(I7,4G15.5,4G15.5)
  250 FORMAT(2X,8('-'),'|',56('-'),'|',56('-'),'|',/,' TOTALS:',
     1 4G15.5,4G15.5)
C
      END
      SUBROUTINE SGWF2SUB7A(HAQ,TLED,CI,SSE,SSV,DZ,DH,DHP,DHC,NN)
C     ******************************************************************
C        ASSEMBLE COEFFICIENTS FOR SOLVING FOR HEAD DISTRIBUTION
C        IN ONE STRING OF CELLS REPRESENTING ONE-HALF OF A DOUBLY
C        DRAINING INTERBED
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GWFSUBMODULE ,ONLY:A1,A2,BB
      DIMENSION DH(NN),DHP(NN),DHC(NN)
C     ------------------------------------------------------------------
C
C1------INITIALIZE
C
      CI2=CI*2.0
      DS=DZ*TLED
      NN1=NN-1
C2------SET COEFFICIENTS FOR CELL BORDERING AQUIFER
      SS=SSE
      A2(1)=CI
      HD=DH(1)
      HC=DHC(1)
      IF(HD.LT.HC) SS=SSV
      A1(1)=-3.*CI-DS*SS
      BB(1)=DS*(SSE*(HC-DHP(1))-HC*SS)-CI2*HAQ-A1(1)*HD
C3------SET COEFFICIENTS FOR INTERIOR CELLS
      BB(2)=-CI*HD
      DO 10 N=2,NN1
      SS=SSE
      A2(N)=CI
      HD=DH(N)
      HC=DHC(N)
      IF(HD.LT.HC) SS=SSV
      CHN=-CI*HD
      BB(N-1)=BB(N-1)+CHN
      BB(N+1)=CHN
      A1(N)=-CI2-DS*SS
      BB(N)=BB(N)+DS*(SSE*(HC-DHP(N))-HC*SS)-A1(N)*HD
   10 CONTINUE
C4------SET COEFFICIENTS FOR CELL BORDERING MIDPLANE OF INTERBED
      SS=SSE
      A2(NN)=CI
      HD=DH(NN)
      HC=DHC(NN)
      BB(NN1)=BB(NN1)-CI*HD
      IF(HD.LT.HC) SS=SSV
      A1(NN)=-CI-0.5*DS*SS
      BB(NN)=BB(NN)+DS*0.5*(SSE*(HC-DHP(NN))-HC*SS)-A1(NN)*HD
C5------RETURN
      RETURN
      END
      SUBROUTINE SGWF2SUB7S(NN)
C     ******************************************************************
C        SOLVE SYSTEM OF EQUATIONS WITH A SYMMETRICAL TRI-DIAGONAL
C        COEFFICIENT MATRIX
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GWFSUBMODULE ,ONLY: A1,A2,BB
C     ------------------------------------------------------------------
C
C1------TRIANGULARIZE LEFT-HAND SIDE MATRIX
      NN1=NN-1
      DO 30 N=1,NN1
      F=1./A1(N)
      C=A2(N)*F
      I=N+1
      A1(I)=A1(I)-C*A2(N)
      A2(N)=C
      BB(I)=BB(I)-C*BB(N)
      BB(N)=BB(N)*F
   30 CONTINUE
      BB(NN)=BB(NN)/A1(NN)
C
C2------BACK SUBSTITE FOR SOLUTION
      DO 40 N=NN1,1,-1
      BB(N)=BB(N)-A2(N)*BB(N+1)
   40 CONTINUE
      RETURN
      END
      SUBROUTINE GWF2SUB72D1D(BUFF,NCOL,NROW,D,ND,LOC)
C     ******************************************************************
C     Move 2-D array into 1-D array
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      DIMENSION BUFF(NCOL,NROW),D(ND)
C     ------------------------------------------------------------------
      L=LOC-1
      DO 10 I=1,NROW
      DO 10 J=1,NCOL
      L=L+1
      D(L)=BUFF(J,I)
   10 CONTINUE
      RETURN
      END
      SUBROUTINE GWF2SUB7SV(IGRID)
C     ******************************************************************
C     SAVE INTERBED STORAGE DATA FOR FUTURE RESTART
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GWFSUBMODULE ,ONLY: DH,DHC,NDF,ND2,IDSAVE
C     ------------------------------------------------------------------
      CALL SGWF2SUB7PNT(IGRID)
C
C1-----PROCESS IF SAVE OPTION SELECTED AND DELAY INTERBEDS EXIST
      IF(IDSAVE.GT.0.AND.NDF) THEN
C2-----WRITE OUT NUMBER OF NODAL HEAD AND PRECONSOLIDATION HEAD VALUES
C2-----THAT ARE BEING SAVED TO DISK
       WRITE(IDSAVE) ND2
C3-----WRITE ARRAYS FOR EACH SYSTEM OF INTERBEDS
       WRITE(IDSAVE) (DH(N),N=1,ND2)
       WRITE(IDSAVE) (DHC(N),N=1,ND2)
      ENDIF
C4-----RETURN
      RETURN
      END

      SUBROUTINE GWF2SUB7DA(IGRID)
C
C     ******************************************************************
C     DEALLOCATE DYNAMIC STORAGE FOR SUB PACKAGE
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GWFSUBMODULE
C     ------------------------------------------------------------------
C
      DEALLOCATE (GWFSUBDAT(IGRID)%IIBSCB)
      DEALLOCATE (GWFSUBDAT(IGRID)%ITMIN)
      DEALLOCATE (GWFSUBDAT(IGRID)%NNDB)
      DEALLOCATE (GWFSUBDAT(IGRID)%NDB)
      DEALLOCATE (GWFSUBDAT(IGRID)%NMZ)
      DEALLOCATE (GWFSUBDAT(IGRID)%NN)
      DEALLOCATE (GWFSUBDAT(IGRID)%ND2)
      DEALLOCATE (GWFSUBDAT(IGRID)%IDSAVE)
      DEALLOCATE (GWFSUBDAT(IGRID)%AC1)
      DEALLOCATE (GWFSUBDAT(IGRID)%AC2)
      DEALLOCATE (GWFSUBDAT(IGRID)%NDF)
      DEALLOCATE (GWFSUBDAT(IGRID)%NNDF)
      DEALLOCATE (GWFSUBDAT(IGRID)%ISBOCF)
      DEALLOCATE (GWFSUBDAT(IGRID)%ISBOCU)
      DEALLOCATE (GWFSUBDAT(IGRID)%OCFLGS)
      DEALLOCATE (GWFSUBDAT(IGRID)%OCLAY)
      DEALLOCATE (GWFSUBDAT(IGRID)%ILSYS)
      DEALLOCATE (GWFSUBDAT(IGRID)%NTSSUM)
      DEALLOCATE (GWFSUBDAT(IGRID)%LN)
      DEALLOCATE (GWFSUBDAT(IGRID)%LDN)
      DEALLOCATE (GWFSUBDAT(IGRID)%NZ)
      DEALLOCATE (GWFSUBDAT(IGRID)%RNB)
      DEALLOCATE (GWFSUBDAT(IGRID)%DH)
      DEALLOCATE (GWFSUBDAT(IGRID)%DHP)
      DEALLOCATE (GWFSUBDAT(IGRID)%DHC)
      DEALLOCATE (GWFSUBDAT(IGRID)%DZ)
      DEALLOCATE (GWFSUBDAT(IGRID)%HC)
      DEALLOCATE (GWFSUBDAT(IGRID)%SCE)
      DEALLOCATE (GWFSUBDAT(IGRID)%SCV)
      DEALLOCATE (GWFSUBDAT(IGRID)%DCOM)
      DEALLOCATE (GWFSUBDAT(IGRID)%A1)
      DEALLOCATE (GWFSUBDAT(IGRID)%A2)
      DEALLOCATE (GWFSUBDAT(IGRID)%BB)
      DEALLOCATE (GWFSUBDAT(IGRID)%SUB)
      DEALLOCATE (GWFSUBDAT(IGRID)%DP)
      DEALLOCATE (GWFSUBDAT(IGRID)%DVB)

C2-----RETURN
      RETURN
      END
      SUBROUTINE SGWF2SUB7PNT(IGRID)
C  Change SUB data to a different grid.
      USE GWFSUBMODULE
C
      IIBSCB=>GWFSUBDAT(IGRID)%IIBSCB
      ITMIN=>GWFSUBDAT(IGRID)%ITMIN
      NNDB=>GWFSUBDAT(IGRID)%NNDB
      NDB=>GWFSUBDAT(IGRID)%NDB
      NMZ=>GWFSUBDAT(IGRID)%NMZ
      NN=>GWFSUBDAT(IGRID)%NN
      ND2=>GWFSUBDAT(IGRID)%ND2
      IDSAVE=>GWFSUBDAT(IGRID)%IDSAVE
      AC1=>GWFSUBDAT(IGRID)%AC1
      AC2=>GWFSUBDAT(IGRID)%AC2
      NDF=>GWFSUBDAT(IGRID)%NDF
      NNDF=>GWFSUBDAT(IGRID)%NNDF
      ISBOCF=>GWFSUBDAT(IGRID)%ISBOCF
      ISBOCU=>GWFSUBDAT(IGRID)%ISBOCU
      OCFLGS=>GWFSUBDAT(IGRID)%OCFLGS
      OCLAY=>GWFSUBDAT(IGRID)%OCLAY
      ILSYS=>GWFSUBDAT(IGRID)%ILSYS
      NTSSUM=>GWFSUBDAT(IGRID)%NTSSUM
      LN=>GWFSUBDAT(IGRID)%LN
      LDN=>GWFSUBDAT(IGRID)%LDN
      NZ=>GWFSUBDAT(IGRID)%NZ
      RNB=>GWFSUBDAT(IGRID)%RNB
      DH=>GWFSUBDAT(IGRID)%DH
      DHP=>GWFSUBDAT(IGRID)%DHP
      DHC=>GWFSUBDAT(IGRID)%DHC
      DZ=>GWFSUBDAT(IGRID)%DZ
      HC=>GWFSUBDAT(IGRID)%HC
      SCE=>GWFSUBDAT(IGRID)%SCE
      SCV=>GWFSUBDAT(IGRID)%SCV
      DCOM=>GWFSUBDAT(IGRID)%DCOM
      A1=>GWFSUBDAT(IGRID)%A1
      A2=>GWFSUBDAT(IGRID)%A2
      BB=>GWFSUBDAT(IGRID)%BB
      SUB=>GWFSUBDAT(IGRID)%SUB
      DP=>GWFSUBDAT(IGRID)%DP
      DVB=>GWFSUBDAT(IGRID)%DVB
C
      RETURN
      END
      SUBROUTINE SGWF2SUB7PSV(IGRID)
C  Save SUB data for a grid.
      USE GWFSUBMODULE
C
      GWFSUBDAT(IGRID)%IIBSCB=>IIBSCB
      GWFSUBDAT(IGRID)%ITMIN=>ITMIN
      GWFSUBDAT(IGRID)%NNDB=>NNDB
      GWFSUBDAT(IGRID)%NDB=>NDB
      GWFSUBDAT(IGRID)%NMZ=>NMZ
      GWFSUBDAT(IGRID)%NN=>NN
      GWFSUBDAT(IGRID)%ND2=>ND2
      GWFSUBDAT(IGRID)%IDSAVE=>IDSAVE
      GWFSUBDAT(IGRID)%AC1=>AC1
      GWFSUBDAT(IGRID)%AC2=>AC2
      GWFSUBDAT(IGRID)%NDF=>NDF
      GWFSUBDAT(IGRID)%NNDF=>NNDF
      GWFSUBDAT(IGRID)%ISBOCF=>ISBOCF
      GWFSUBDAT(IGRID)%ISBOCU=>ISBOCU
      GWFSUBDAT(IGRID)%OCFLGS=>OCFLGS
      GWFSUBDAT(IGRID)%OCLAY=>OCLAY
      GWFSUBDAT(IGRID)%ILSYS=>ILSYS
      GWFSUBDAT(IGRID)%NTSSUM=>NTSSUM
      GWFSUBDAT(IGRID)%LN=>LN
      GWFSUBDAT(IGRID)%LDN=>LDN
      GWFSUBDAT(IGRID)%NZ=>NZ
      GWFSUBDAT(IGRID)%RNB=>RNB
      GWFSUBDAT(IGRID)%DH=>DH
      GWFSUBDAT(IGRID)%DHP=>DHP
      GWFSUBDAT(IGRID)%DHC=>DHC
      GWFSUBDAT(IGRID)%DZ=>DZ
      GWFSUBDAT(IGRID)%HC=>HC
      GWFSUBDAT(IGRID)%SCE=>SCE
      GWFSUBDAT(IGRID)%SCV=>SCV
      GWFSUBDAT(IGRID)%DCOM=>DCOM
      GWFSUBDAT(IGRID)%A1=>A1
      GWFSUBDAT(IGRID)%A2=>A2
      GWFSUBDAT(IGRID)%BB=>BB
      GWFSUBDAT(IGRID)%SUB=>SUB
      GWFSUBDAT(IGRID)%DP=>DP
      GWFSUBDAT(IGRID)%DVB=>DVB
C
      RETURN
      END
